#Andrew Gooch
#April 11, 2017
#Ripping Yarn, Political Communications
#Replication R file 1: Making Figure 1

#This R file created two panels using ggplot2 and then merges them together to make one figure

#required packages
library(MASS)
library(arm)
library(ggplot2)
library(grid)
library(gridExtra)
library(foreign)
library(survey)
library(plotrix)

#loads data. Please change to your file directory where the data is located
Fig1a <- read.csv("C:/Users/agooch/Dropbox/Working papers/Ripping Yarn/Data/Convention count_all_appeals_aggregated_by_year_FIG1a.csv")
attach(Fig1a)
head(Fig1a)

#Percent Storytelling out of all appeals
n = all_appeals    # valid responses count 
k = storytelling
avg_st = k/n
avg_st #average storytelling per year as a vector
mean(avg_st) #avg storytelling from 1964-2012 (25.8%)


str(prop.test(k[year==2012], n[year==2012]) )

#use prop test function to get upper and lower CI bounds
#lower bounds by year
lb_2012 = prop.test(k[year==2012], n[year==2012])$conf.int[1:1]
lb_2008 = prop.test(k[year==2008], n[year==2008])$conf.int[1:1]
lb_2004 = prop.test(k[year==2004], n[year==2004])$conf.int[1:1]
lb_2000 = prop.test(k[year==2000], n[year==2000])$conf.int[1:1]
lb_1996 = prop.test(k[year==1996], n[year==1996])$conf.int[1:1]
lb_1992 = prop.test(k[year==1992], n[year==1992])$conf.int[1:1]
lb_1988 = prop.test(k[year==1988], n[year==1988])$conf.int[1:1]
lb_1984 = prop.test(k[year==1984], n[year==1984])$conf.int[1:1]
lb_1980 = prop.test(k[year==1980], n[year==1980])$conf.int[1:1]
lb_1976 = prop.test(k[year==1976], n[year==1976])$conf.int[1:1]
lb_1972 = prop.test(k[year==1972], n[year==1972])$conf.int[1:1]
lb_1968 = prop.test(k[year==1968], n[year==1968])$conf.int[1:1]
lb_1964 = prop.test(k[year==1964], n[year==1964])$conf.int[1:1]

#upper bounds by year
ub_2012 = prop.test(k[year==2012], n[year==2012])$conf.int[2:2]
ub_2008 = prop.test(k[year==2008], n[year==2008])$conf.int[2:2]
ub_2004 = prop.test(k[year==2004], n[year==2004])$conf.int[2:2]
ub_2000 = prop.test(k[year==2000], n[year==2000])$conf.int[2:2]
ub_1996 = prop.test(k[year==1996], n[year==1996])$conf.int[2:2]
ub_1992 = prop.test(k[year==1992], n[year==1992])$conf.int[2:2]
ub_1988 = prop.test(k[year==1988], n[year==1988])$conf.int[2:2]
ub_1984 = prop.test(k[year==1984], n[year==1984])$conf.int[2:2]
ub_1980 = prop.test(k[year==1980], n[year==1980])$conf.int[2:2]
ub_1976 = prop.test(k[year==1976], n[year==1976])$conf.int[2:2]
ub_1972 = prop.test(k[year==1972], n[year==1972])$conf.int[2:2]
ub_1968 = prop.test(k[year==1968], n[year==1968])$conf.int[2:2]
ub_1964 = prop.test(k[year==1964], n[year==1964])$conf.int[2:2]

#pull each year average from the vector
m_2012 = avg_st[year==2012]
m_2008 = avg_st[year==2008]
m_2004 = avg_st[year==2004]
m_2000 = avg_st[year==2000]
m_1996 = avg_st[year==1996]
m_1992 = avg_st[year==1992]
m_1988 = avg_st[year==1988] 
m_1984 = avg_st[year==1984]
m_1980 = avg_st[year==1980]
m_1976 = avg_st[year==1976]
m_1972 = avg_st[year==1972]
m_1968 = avg_st[year==1968]
m_1964 = avg_st[year==1964]

lb <-c(lb_1964*100,
       lb_1968*100,lb_1972*100,
       lb_1976*100,lb_1980*100,
       lb_1984*100,lb_1988*100,
       lb_1992*100,lb_1996*100,
       lb_2000*100,lb_2004*100,
       lb_2008*100,lb_2012*100)

ub <-c(ub_1964*100,
       ub_1968*100,ub_1972*100,
       ub_1976*100,ub_1980*100,
       ub_1984*100,ub_1988*100,
       ub_1992*100,ub_1996*100,
       ub_2000*100,ub_2004*100,
       ub_2008*100,ub_2012*100)


coef <- c(m_1964*100,
          m_1968*100,m_1972*100,
          m_1976*100,m_1980*100,
          m_1984*100,m_1988*100,
          m_1992*100,m_1996*100,
          m_2000*100,m_2004*100,
          m_2008*100,m_2012*100)

yearly_num <- c(1,2,3,4,5,6,7,8,9,10,11,12,13)

#calculate a fitted line of storytelling over time

lm(coef ~ yearly_num) # intercept = 17.922 and slope = 1.127




#turns vectors into dataframe for ggplot2
coef_df<-data.frame(coef) 
row.names(coef_df) <- c("1964","1968", "1972", "1976", "1980", "1984", "1988", "1992", "1996", "2000", "2004", "2008", "2012")


df<-data.frame(Intercepts=coef,lb,ub,
               lev.names= rownames(coef_df))


p1 <- ggplot(df,aes(lev.names,Intercepts,shape=lev.names))+ 
  geom_errorbar(aes(ymin=lb, ymax=ub), width=0,color="gray10", lty=1, size=1) + 
  geom_point(aes(size=2)) +
  geom_abline(intercept = 17.922, slope = 1.127) +
  scale_color_manual(name=" ",values=c("black","darkgrey")) + 
  scale_shape_manual(name=" ",values=c(16,16,16,16,16,16,
                                       16,16,16,16,16,16,
                                       16)) + 
  theme_bw() + 
  xlab(" ") + 
  ylab("%") +
  ggtitle("Percentage of Storytelling Appeals") +
  guides(size=FALSE,shape=FALSE) +
  theme(axis.text.x=element_text(size=rel(1.6)),
        axis.title.x=element_text(size=rel(1.6)),
        axis.title.y=element_text(size=rel(1.4)),
        axis.text.y=element_text(size=rel(1.6)),
        panel.grid.minor=element_blank(),
        panel.grid.major.x=element_blank(),
        axis.ticks.y = element_blank(), 
        axis.ticks.x = element_blank(), 
        panel.border = element_blank(),
        axis.line.x = element_line(color = 'black') )


print(p1) #top panel of FIG 1


################################################



Fig1b <- read.csv("C:/Users/agooch/Dropbox/Working papers/Ripping Yarn/Data/Storytelling_breakdown_aggregated_by_year_final_FIG1b.csv")
attach(Fig1b)
head(Fig1b)

p2 <- ggplot(Fig1b, aes(YEAR)) + 
  geom_line(aes(y = IMS_per, colour = "Impersonalized"), linetype = 2, size=1.25) + 
  geom_line(aes(y = BIO_per, colour = "Autobiographical"), linetype = 3, size=1.25)+
  geom_line(aes(y = PS_per, colour = "Personalized"), linetype = 1, size=1.25)+
  scale_x_continuous(breaks=seq(1964, 2012, 4))+
  theme_bw() + 
  xlab(" ") + 
  ylab("%") +
  ggtitle("Storytelling Appeals Only (%)") +
  guides(size=FALSE,shape=FALSE) +
  theme(axis.text.x=element_text(size=rel(1.6)),
        axis.title.x=element_text(size=rel(1.6)),
        axis.title.y=element_text(size=rel(1.4)),
        axis.text.y=element_text(size=rel(1.6)),
        panel.grid.minor=element_blank(),
        panel.grid.major.x=element_blank(),
        axis.ticks.y = element_blank(), 
        axis.ticks.x = element_blank(), 
        panel.border = element_blank(),
        axis.line.x = element_line(color = 'black'),
        legend.position = c(.75, 0.75),
        legend.title=element_blank(),
        legend.text = element_text(size=rel(1.4)))

print(p2)



p3 <- grid.arrange(arrangeGrob(p1,p2,nrow=2)) 






